Capítulo 5 Analise geoespacial
5.1 Dados espaciais e geoespaciais
Dados espaciais são os que utilizam o sistema de coordenadas cartesianas com três (x, y e z) ou mais dimensões. Dados geoespaciais são dados que podem ser mapeados no planeta Terra e relacionadas com outros dados baseados em sistemas de coordenadas geográficas. Como grande parte dos dados ambientais são afetados pela localização geográfica, a análise geoespacial traz informações importantes.
A teoria por trás da análise feoespacial já é coberta pela disciplina de Cartografia e Geoprocessamento, portanto, aqui iremos abordar somente a parte prática.
No capítulo, serão utilizados os seguintes pacotes: sf para ler e trabalhar com dados espaciais e mapview para a etapa de criação de mapas.
No curso de Engenharia Ambiental e Urbana, utilizam-se shapefiles (.shp) para realizar análises espaciais, portanto, estes serão utilizados no capítulo.
5.1.1 Sistema Geodesico de Referência (SGR)
Basicamente, o SGR é um sistema de coordenadas definido a partir de um elipsóide de referência, posicionado e orientado em relação à superfície da Terra. A partir dele, é possível localizar espacialmente qualquer feição na superfície terrestre. Os mais conhecidos são: SAD69, WGS84 e o SIRGAS 2000.
Aplicação
Para a aplicação será reproduzido o exemplo do TBEP R Training. Para isso, serão instalados os pacotes sf e mapview.
options(repos = list(CRAN="http://cran.rstudio.com/"))
options("install.lock"=FALSE)
install.packages(c('sf','mapview'))## Installing packages into 'C:/Users/Luiz Arthur/Dropbox/PC/Documents/R/win-library/4.1'
## (as 'lib' is unspecified)
## package 'sf' successfully unpacked and MD5 sums checked
## package 'mapview' successfully unpacked and MD5 sums checked
##
## The downloaded binary packages are in
## C:\Users\Luiz Arthur\AppData\Local\Temp\RtmpCcWuwO\downloaded_packages
library(sf)## Warning: package 'sf' was built under R version 4.1.3
## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE
library(mapview)## Warning: package 'mapview' was built under R version 4.1.3
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5 v purrr 0.3.4
## v tibble 3.1.6 v dplyr 1.0.8
## v tidyr 1.2.0 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1
## Warning: package 'dplyr' was built under R version 4.1.3
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
O shapefile “sgdat.shp” são dados da cobertura de algas marinhas em Tampa Bay em 2016. As “features” são as linhas do vetor e os “fields” são as colunas, ou melhor, atributos (“OBJECT ID” e “FLUCCS”). O SGR do arquivo é WGS 84. A coluna “geometry” armazena os dados espaciais (longitude e latitude).
#sgdat shapefile
sgdat <- st_read('Data/sgdat.shp')## Reading layer `sgdat' from data source
## `C:\Users\Luiz Arthur\Dropbox\PC\Documents\UFABC Beatriz\TG Beatriz Lima\TG_Beatriz_Lima\Data\sgdat.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 575 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -82.72462 ymin: 27.81386 xmax: -82.48426 ymax: 28.03548
## Geodetic CRS: WGS 84
# utilidades do pacote sf
methods(class="sf")## [1] $<- [
## [3] [[<- aggregate
## [5] anti_join arrange
## [7] as.data.frame cbind
## [9] coerce dbDataType
## [11] dbWriteTable distinct
## [13] dplyr_reconstruct filter
## [15] full_join gather
## [17] group_by group_split
## [19] identify initialize
## [21] inner_join left_join
## [23] mapView merge
## [25] mutate nest
## [27] pivot_longer pivot_wider
## [29] plot print
## [31] rbind rename
## [33] right_join rowwise
## [35] sample_frac sample_n
## [37] select semi_join
## [39] separate separate_rows
## [41] show slice
## [43] slotsFromS3 spread
## [45] st_agr st_agr<-
## [47] st_area st_as_s2
## [49] st_as_sf st_as_sfc
## [51] st_bbox st_boundary
## [53] st_buffer st_cast
## [55] st_centroid st_collection_extract
## [57] st_convex_hull st_coordinates
## [59] st_crop st_crs
## [61] st_crs<- st_difference
## [63] st_drop_geometry st_filter
## [65] st_geometry st_geometry<-
## [67] st_inscribed_circle st_interpolate_aw
## [69] st_intersection st_intersects
## [71] st_is st_is_valid
## [73] st_join st_line_merge
## [75] st_m_range st_make_valid
## [77] st_minimum_rotated_rectangle st_nearest_points
## [79] st_node st_normalize
## [81] st_point_on_surface st_polygonize
## [83] st_precision st_reverse
## [85] st_sample st_segmentize
## [87] st_set_precision st_shift_longitude
## [89] st_simplify st_snap
## [91] st_sym_difference st_transform
## [93] st_triangulate st_union
## [95] st_voronoi st_wrap_dateline
## [97] st_write st_z_range
## [99] st_zm summarise
## [101] transform transmute
## [103] ungroup unite
## [105] unnest
## see '?methods' for accessing help and source code
Esse é o passo a passo de como importar um shapefile. Porém, muitas vezes não possuímos um shapefile e queremos criar um a partir de um dataframe. Para isso, é necessário que o dataframe inclua as coordenadas geográficas (longitude e latitude) e que tenhamos conhecimento do SGR. O dataframe ´fishdat´ possui as características dos peixes encontrados e o statloc apresenta a localização deles. O passo a passo será realizado abaixo.
getwd()## [1] "C:/Users/Luiz Arthur/Dropbox/PC/Documents/UFABC Beatriz/TG Beatriz Lima/TG_Beatriz_Lima"
# dados da presença de peixes em Tampa Bay
fishdat <- read.csv("Data/fishdat.csv")
statloc <- read.csv("Data/statloc.csv")# estrutura dos dados
str(fishdat)## 'data.frame': 2844 obs. of 12 variables:
## $ OBJECTID : int 1550020 1550749 1550750 1550762 1550828 1550838 1550842 1551131 1551311 1551335 ...
## $ Reference : chr "TBM1996032006" "TBM1996032004" "TBM1996032004" "TBM1996032207" ...
## $ Sampling_Date: chr "1996-03-20" "1996-03-20" "1996-03-20" "1996-03-22" ...
## $ yr : int 1996 1996 1996 1996 1996 1996 1996 1996 1996 1996 ...
## $ Gear : int 300 22 22 20 160 300 300 300 300 22 ...
## $ ExDate : chr "2018-04-12 10:27:38" "2018-04-12 10:25:23" "2018-04-12 10:25:23" "2018-04-12 10:25:23" ...
## $ Bluefish : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Common.Snook : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Mullets : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Pinfish : int 0 54 0 80 0 0 0 0 1 1 ...
## $ Red.Drum : int 0 0 1 0 4 0 0 0 0 0 ...
## $ Sand.Seatrout: int 1 0 0 0 0 1 5 66 0 0 ...
str(statloc)## 'data.frame': 2173 obs. of 3 variables:
## $ Reference: chr "TBM1996032006" "TBM1996032004" "TBM1996032207" "TBM1996042601" ...
## $ Latitude : num 27.9 27.9 27.9 28 27.9 ...
## $ Longitude: num -82.6 -82.6 -82.5 -82.7 -82.6 ...
Para isso, utilizaremos a função st_as_sf() para transformar o dataframe em um objeto sf. Primeiramente, precisamos juntar os dois datasets (fishdat e statloc) e dizer qual coluna que possui os dados da geometria (latitude e longitude). Além disso, é necessário dizer qual o SGR e, além disso, precisamos garantir que ambos datasets possuam o mesmo SGR. Por enquanto, podemos fazer um “chute calibrado” que é o WGS84.
#juntando os dois dataframes
joindata <- left_join(fishdat,statloc,by="Reference")
#criando o objeto de dados espaciais
joindata <- st_as_sf(joindata, coords=c('Longitude','Latitude'), crs = st_crs(sgdat))
#tipo de objeto sf
str(joindata)## Classes 'sf' and 'data.frame': 2844 obs. of 13 variables:
## $ OBJECTID : int 1550020 1550749 1550750 1550762 1550828 1550838 1550842 1551131 1551311 1551335 ...
## $ Reference : chr "TBM1996032006" "TBM1996032004" "TBM1996032004" "TBM1996032207" ...
## $ Sampling_Date: chr "1996-03-20" "1996-03-20" "1996-03-20" "1996-03-22" ...
## $ yr : int 1996 1996 1996 1996 1996 1996 1996 1996 1996 1996 ...
## $ Gear : int 300 22 22 20 160 300 300 300 300 22 ...
## $ ExDate : chr "2018-04-12 10:27:38" "2018-04-12 10:25:23" "2018-04-12 10:25:23" "2018-04-12 10:25:23" ...
## $ Bluefish : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Common.Snook : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Mullets : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Pinfish : int 0 54 0 80 0 0 0 0 1 1 ...
## $ Red.Drum : int 0 0 1 0 4 0 0 0 0 0 ...
## $ Sand.Seatrout: int 1 0 0 0 0 1 5 66 0 0 ...
## $ geometry :sfc_POINT of length 2844; first list element: 'XY' num -82.6 27.9
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:12] "OBJECTID" "Reference" "Sampling_Date" "yr" ...
#checando SGR
st_crs(joindata)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
st_crs(sgdat) ## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
Caso seja necessário modificar a projeção, utiliza-se a função `st_transform(). Nesse caso, não precisamos modificar já que o shapefile (sgdat) tem o mesmo SGR do que estamo querendo criar.
Agora, iniciaremos a análise geoespacial dos dados. Inicialmente, iremos dar uma olhada geral para entender qual os dados que estamos lidando.
O padrão é que a função ´plot()´ plote todas as feições. Para plotar somente a geometria, utiliza-se st_geometry().
plot(st_geometry(joindata)) 
plot(joindata)## Warning: plotting the first 10 out of 12 attributes; use max.plot = 12 to plot
## all

plot(sgdat)
Conforme observamos o shapefile “sgdat” com os dados das algas marinhas e o “joindata” com os dados do posicionamento de peixes, é possível verificar que existem áreas de intersecção entre ambos. Para analisar novamente, iremos plotar somente a geometria de ambos:
plot(joindata$geometry)
plot(sgdat$geometry)
Vamos filtrar somente os dados dos peixes do ano de 2016:
filt_data <- joindata %>%
filter(yr == 2016)
plot(st_geometry(filt_data))
Agora, verificaremos quantos peixes foram vistos nos mesmos locais em que encontraram-se algas marinhas em 2016. Ou seja, iremos selecionar as localizações que possuem ambos dados. Para isso, iremos utilizar o código abaixo:
fish_crop <- filt_data[sgdat, ]
plot(fish_crop$geometry)
O que foi realizado até agora é somente a intersecção da geometria de ambos datasets. Portanto, agora realizaremos a intersecção de ambos dados, incluindo atributos:
fish_int <- st_intersection(filt_data, sgdat)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(st_geometry(fish_int))
View(fish_int)É possível utilizar ferramentas do tidyverse. Abaixo, iremos fazer a soma de todos os Pinfish foram pegos em 2016:
fish_cnt <- fish_int %>%
group_by(FLUCCS) %>%
summarise(
cnt = sum(Pinfish)
)
fish_cnt## Simple feature collection with 2 features and 2 fields
## Geometry type: MULTIPOINT
## Dimension: XY
## Bounding box: xmin: -82.7182 ymin: 27.82623 xmax: -82.53237 ymax: 28.02418
## Geodetic CRS: WGS 84
## # A tibble: 2 x 3
## FLUCCS cnt geometry
## <chr> <int> <MULTIPOINT [°]>
## 1 9113 1559 ((-82.53352 27.92907), (-82.53617 27.91043), (-82.53535 27.91212~
## 2 9116 4766 ((-82.55893 27.96612), (-82.5588 27.96633), (-82.55797 27.96632)~
Além de realizar a soma nos atributos numéricos (quantidade de Pinfish), também é realizada nos atributos geométricos (latitude e longitude). Conforme apresentado na tabela anterior, existe uma maior quantidade de Pinfishs em áreas onde existe maior quantidade de algas marinhas (FLUCCS=9116). É possível realizar um gráfico em relação às duas categorias de cobertura de algas marinhas (´9113´: desigual, ´9116´: contínua).
ggplot(fish_cnt, aes(x = FLUCCS, y = cnt)) +
geom_bar(stat = 'identity', fill='navyblue')
Agora será realizada a confecção de mapas. Utilizaremos os pacotes ggplot2 inicialmente:
ggplot() +
geom_sf(data = sgdat, fill = 'green') +
geom_sf(data = joindata)
Agora, para criar um mapa interativo para selecionar e dar zoom nos dados, utilizaremos o pacote mapview:
mapview(sgdat, col.regions = 'green') +
mapview(joindata, zcol = 'Gear')